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We study attractively interacting spin-| fermions on the square lattice subject to a spin population 
imbalance. Using unbiased diagrammatic Monte Carlo simulations we find an extended region in 
the parameter space where the Fermi liquid is unstable towards formation of Cooper pairs with non¬ 
zero center-of-mass momentum, known as the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state. In 
contrast to earlier mean-field and quasi-classical studies we provide quantitative and well-controlled 
predictions on the existence and location of the relevant Fermi-liquid instabilities. The highest 
temperature where the FFLO instability can be observed is about half of the superfluid transition 
temperature in the unpolarized system. 
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Fifty years after the initial prediction by Fnlde, Ferrell, 
Larkin, and Ovchinnikov (FFLO) [Il[2], superconducting 
phases with spontaneously broken translational invari¬ 
ance are still at the center of interest in such diverse fields 
as solid state physics, cold atomic gases, nuclear physics, 
and dense quark matter in neutron stars [3H7]. While 
the underlying mechanism is generic enough to apply to 
any partially polarized Fermi system, it has proven sur¬ 
prisingly difficult to unambiguously observe such phases 
in nature. Recently, however, experimental evidence has 
been mounting for their existence in heavy fermion com¬ 
pounds HHini and layered organic materials pThT^ . On 
the other hand, experiments with ultracold atoms, which 
are among the cleanest imbalanced Fermi systems with¬ 
out the need for a magnetic field, so far failed to demon¬ 
strate inhomogeneous superfluidity [igin] — although 
there is some evidence for such a phase in one dimension 
(ID) m — possibly due to small extent of the param¬ 
eter region where an FFLO phase may exist in three di¬ 
mensions (3D) and difficulties in reaching sufficiently low 
temperatures [5]. 

On the theoretical side, results on the existence and 
nature of FFLO phases based on well-controlled micro¬ 
scopic theories are scarce, with the exception of ID sys¬ 
tems, where exact analytical and numerical studies are 
possible [T9]-[25], and where finite-momentum pairing is 
a generic feature of the spin-imbalanced phase diagram. 
In higher dimensions, most studies are based on effective 
field theories in the neighborhood of critical points or re¬ 
sort to quasi-classical or mean-field approximations. For 
3D Fermi gases, many features of the mean-field phase 
diagram [26] have been corroborated by fixed-node diffu¬ 
sion quantum Monte Carlo calculations m-, whether the 
FFLO phase does exist in a small sliver of the phase di¬ 
agram, as predicted by the mean-field theory, is however 
still subject to debate. 

The FFLO state is expected EE] to occupy a larger 


parameter region in two dimensions (2D), and lattice 
effects may further increase its stability [28l [29]. Cor¬ 
respondingly, mean-field calculations |30l [31] and real- 
space dynamical mean-field theory (DMFT) for fermions 
in anisotropic optical lattices find a stable and extended 
spatially modulated superfluid [321(34] . However, such 
approximations are particularly questionable in 2D. The 
only numerically exact study to date is a determinan- 
tal quantum Monte Carlo simulation of the attractive 
Hubbard model [35], showing a finite-momentum peak 
in the pair-momentum distribution in large parts of the 
polarization-temperature phase diagram. Unfortunately, 
this study is severely limited by the negative sign problem 
and could not reach low enough temperatures to establish 
phase coherence of the pairs. Therefore, the question of 
whether or not an FFLO phase with (quasi-)long-range 
order can emerge in a given microscopic model remains 
open. 

In this Letter we employ the unbiased diagrammatic 
Monte Carlo (DiagMC) method [36H38] to identify su¬ 
perfluid instabilities of spin-imbalanced fermions on a 2D 
lattice in a controlled way at lower temperatures than 
hitherto accessible. Our main result is that, for attractive 
interactions of the order of half the bandwidth, there is an 
extended region in the temperature-polarization plane, 
indicated by red shading in Fig. [^ where the Fermi liq¬ 
uid is unstable exclusively towards FFLO superfluidity. 
Specifically, we simulate the Hubbard model on a square 
lattice 

H = -t + h.c^ +uY^ (1) 

i 

with nearest-neighbor hopping amplitude t = 1 setting 
the scale of energy and on-site attraction U < 0. Here 
and c-^ are fermionic creation and annihilation operators 
with spin a =t, and riia- = Spin imbalance is 

quantified by the polarization P = 


2 



FIG. 1. (all figures: color online). Phase diagram for 
Ujt — —4 at quarter filling: The white region is a Fermi 
liquid. In the blue shaded region, the Fermi liquid is unsta¬ 
ble towards conventional (Q = 0) pairing. In the red shaded 
region there is an exclusive FFLO instability with finite pair 
momentum Q^. Open symbols indicate whether zero- (blue 
circles) or finite-momentum pairing (red diamonds) is dom¬ 
inant (black squares: no significant difference). The black 
dotted line separates the two regimes. All lines are guides to 
the eye. 


such that P = 1 corresponds to a fully polarized system. 
In the following, we present results for U = —4 at quarter 
filling n = = 0.5. 


Our DiagMC algorithm nil ED stochastically sam¬ 
ples many-body Feynman diagrams (built on the bare 
Green’s function) for the self-energy and the two-particle- 
irreducible pairing vertex directly in the thermodynamic 
limit. Due to the diagrammatic sign problem, in practice 
a cutoff on the maximum addressed diagram order 
is introduced and independence of the results on is 
checked by varying the cutoff. We identify continuous 
phase transitions to ordered phases by monitoring the 
divergence of the corresponding susceptibilities on ap¬ 
proach to the phase boundary from the normal phase. 
According to the Bethe-Salpeter equation, the suscep- 


tibility Xc = cXc = xi°V 




m a 


given channel c (with, e.g., zero or non-zero center-of- 
mass momentum) diverges when the largest eigenvalue 
of the kernel reaches unity. Here, Xc^^ denotes the 

product of two one-particle propagators and Fc the two- 
particle-irreducible pairing vertex, see Refs. [39l [40] for 
details. In the present case, we are primarily interested 
in pairing of t- with ^-particles, which gives rise to both 
the conventional BCS and the FFLO phases, and hence 
concentrate on the superconducting channels with total 
spin projection 5';^ = 0 first. We refer to these as “sin¬ 
glet” channels im and compute their pairing eigenvalues 
Aq for different pair momenta Q. 



FIG. 2. (top) Temperature dependence of the leading pairing 
eigenvalues for zero momentum (open symbols/dashed lines) 
and finite momentum (filled symbols/solid lines) for dif¬ 
ferent polarizations P. (center) Estimates of the critical po¬ 
larization for the two channels from varying diagram order 
cutoff W for temperatures T = 0.5 (black), T = 0.1 (green), 
and T = 0.025 (cyan). The left-most data points show our 
extrapolations W ^ oo. (bottom) Difference between FFLO 
and conventional pairing eigenvalues for varying polarization. 
The corresponding temperatures are indicated to the right of 
the curves. Gircles and diamonds on the curves indicate the 
critical polarization where zero- and finite-momentum eigen¬ 
values respectively cross unity. Also shown is the pair mo¬ 
mentum magnitude Q* (red dots). 


It is possible that the transition to an FFLO state is 
actually first-order due to appearance of solid-type or¬ 
der. In this case, the transition temperature extracted 
from the Bethe-Salpeter equation would correspond to a 
lower bound. However, the FFLO transition on the 2D 
lattice is generally believed to be continuous, at least in 
the neighborhood of the temperature where the FFLO 
instability first emerges [3Q1[3T1[424[44] . 

In the following we compare the pairing eigenvalues at 
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(3 = 0 and at a non-zero candidate momentum (to be 
defined below). Studying their temperature dependence, 
shown in the top panel of Fig. we find that a non-zero 
polarization strongly suppresses the singlet superfluid in¬ 
stabilities as soon as the temperature is low enough to 
resolve the mismatch between the Fermi surfaces (FSs) 
of the minority and majority species: While the tran¬ 
sition temperature in the unpolarized system is roughly 
Tc = 0.15, a moderate polarization of P = 0.2 may only 
lead to a transition (in the FFLO channel) at the lower 
end of the considered temperature range, Tc ^ 0.025. At 
larger polarizations P > 0.3 all the singlet eigenvalues 
seem to saturate below unity, indicating the absence of a 
transition in the considered channels at any temperature. 

Comparing eigenvalues for zero and non-zero pair mo¬ 
mentum, one may differentiate three regimes: At high 
temperatures the FSs are so blurred that the two chan¬ 
nels are essentially degenerate. In the region where the 
effects of the FS mismatch first become noticeable, there 
is a small advantage of zero-momentum pairing. Here 
a configuration where all parts of one FS are close to 
the other FS, even if the two never intersect, is ap¬ 
parently more favorable than the alternative with some 
matching parts and others that are very far apart. At 
an even lower temperature, finally, the zero-momentum 
eigenvalue starts decreasing whereas the finite momen¬ 
tum one continues growing, although with a decreasing 
rate. Depending on polarization, the following scenarios 
are generically realised when the system is cooled down: 
(a) For small polarization, the Q = 0 eigenvalue grows 
to unity before it is overtaken by the eigenvalue, (b) 
For larger polarization, the FFLO eigenvalue will reach 
unity first, (c) For even larger polarization, all singlet 
eigenvalues saturate below unity, i.e. the Fermi liquid 
phase remains stable until triplet pairing emerges at ex¬ 
ponentially low temperatures (see below). 

By tracking the growth of the pairing eigenvalues with 
decreasing polarization at fixed temperature T, we find 
the critical polarization Pc{T) for superfluidity. Its ex¬ 
trapolation in the diagram-order cutoff, the uncertainty 
of which is indicated by horizontal error bars on the 
phase boundaries of Fig. is shown in the central panel 
of Fig. The bottom panel of Fig. shows the dif¬ 
ference between non-zero- and zero-momentum pairing 
eigenvalues. A positive (negative) difference corresponds 
to dominant FFLO (conventional) pairing fluctuations in 
the Fermi liquid phase, indicated by red diamonds (blue 
circles) in the phase diagram. Non-zero-momentum pair¬ 
ing fluctuations are dominant at large polarization and 
low temperature, which is in accord with the large region 
found in Ref. [35] where the pair momentum distribution 
function is peaked at finite momenta. For temperatures 
T < 0.05, the difference is positive at the critical po¬ 
larization Pc{T) implying that the FFLO instability is 
reached before the conventional superfluid one in this 
region of the phase diagram m- We cannot reliably 
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FIG. 3. Influence of density shown for n = 0.8 (top row) and 
n = 0.1 (center row) with polarization P — 0.3 and temper¬ 
ature T = 0.025. Left panels show majority spin momentum 
distribution (colors, from blue=unoccupied to red=occupied) 
and minority FS (dashed contour), as well as the latter shifted 
by the optimal pair momentum (solid contour). The other 
panels illustrate the dependence of the one-particle propaga¬ 
tor product x{Q) on the pair momentum Q. (bottom) De¬ 
pendence of the optimal pair momentum (3*, extracted from 
the one-particle propagator product x((3), on density n and 
polarization P. Dotted lines indicate the weak-coupling form 
for an isotropic dispersion, dashed lines for a square-shaped 
Fermi surface. 


compute the extent of the FFLO phase because our dia¬ 
grammatic approach is not valid inside an ordered phase, 
but we can estimate the region of onset of conventional 
order in the absence of an FFLO phase by extrapolating 
the growth of the corresponding pairing eigenvalue from 
the Fermi liquid phase. As means of extrapolation we use 
the finite-order eigenvalue estimates inside the superfluid 
phase, drawn in lighter shades in Fig. 

In our study, the optimal pair momentum of the FFLO 
state could only be determined approximately be¬ 
cause an optimization of the pairing eigenvalue Aq over 
Q would be too costly within DiagMC. To this end we 
replace the irreducible vertex in the Bethe-Salpeter equa¬ 
tion by the bare interaction I/, such that the approximate 
pairing eigenvalue Aq = —Ux{Q) only depends on Q via 
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the product of single-particle propagators 

X(Q)=r^ / —Gl{k,iLOn)G^{Q-k,-iuJn), (2) 

n 

which can easily be evaluated numerically. The optimal 
pair momentum is always found to lie on the coordinate 
axes of the Brillouin zone. This is most easily under¬ 
stood close to half filling (top row of Fig. |^, where the 
FSs are almost squares. Then, a pair momentum of the 
form = ((5^,0) (and those related by point-group 
symmetry) can connect two sides of the minority FS to 
the corresponding majority FS patches, whereas, say, a 
diagonal pair momentum could only connect one side of 
each FS. For dilute systems (center row of Fig. [^, the 
FSs are almost isotropic such that the majority and mi¬ 
nority FS can at best touch in one tangential point. The 
difference between pair momenta with the same magni¬ 
tude is rather small, but for finite filling we always find 
a slight preference for pair momenta along the lattice 
axes. The bottom panel of Fig. plots this pair mo¬ 
mentum Q* found by numerical optimization for differ¬ 
ent site fillings and polarizations. In general, there is 
no closed expression for but one can consider two 
limiting cases: (a) For circular FSs the respective Fermi 
momenta are kp = so the t and ^ FSs are con¬ 
nected by = \/2^(\/r^FP—\/r^-?). (b) 

For square-shaped FSs, whose corners lie on the coordi¬ 
nate axes at kp = \/27r^^, the optimal pair momentum 
is Q* = \/7r^n(\/l + P — \/l ~ P)- These estimates are 
indicated by dotted and dashed lines, respectively, in the 
bottom panel of Fig. The data obtained by numer¬ 
ical optimization lie quite consistently between the two 
extreme estimates. While the approximation Aq will in 
general strongly overestimate the pairing eigenvalue due 
to the neglect of correlation effects in the vertex, the ex¬ 
tracted pair momentum is expected to be accurate 
because the momentum dependence of the vertex is typi¬ 
cally much weaker than that of the propagators. Strictly 
speaking, we cannot exclude the (unlikely) existence of a 
stronger instability at a different momentum. This means 
that our phase diagram is conservative in the sense that 
the region of the FFLO phase might become only larger 
if additional pair momenta are relevant. 

At strong polarization and at weaker interaction, all 
singlet-pairing eigenvalues saturate below unity. Here, 
triplet pairing, which is not susceptible to the FS mis¬ 
match, will emerge at (exponentially) low temperatures 
[46] due to an effective interaction between identical par¬ 
ticles mediated by the other species, just as in the case 
of a spin-dependent hopping anisotropy [39]. In prin¬ 
ciple, either the majority or the minority species may 
have the dominant instability. In second-order pertur¬ 
bation theory at quarter filling, the majority species al¬ 
ways reach the superfluid transition first, independent 
of the polarization. We have confirmed this with Di- 



FIG. 4. Pairing eigenvalues for different channels at quarter 
filling and polarization P = 0.2. Shown are the singlet pair¬ 
ing eigenvalues with pair momenta Q = 0, (i.e. on the 

coordinate axis) and (o^ diagonal) as well as the 

p-wave triplet eigenvalues for zero-momentum pairing of the 
majority (p^^) and minority (p^^) species, respectively. The 
latter have been multiplied by a factor of ten. 


agMC calculations for P = 0.2 (Fig. and P = 0.4 (not 
shown). In Fig.we compare the eigenvalues in five dif¬ 
ferent channels: Among the singlet-pairing eigenvalues, 
the pair momentum dominates at low temperatures, 
whereas the conventional Q = 0 channel saturates below 
T < 0.2. A further candidate, which is the optimal 
pair momentum on the BZ diagonal within the approx¬ 
imation is always subdominant to Q^. The triplet 
eigenvalues, on the other hand, are by an order of magni¬ 
tude smaller at the temperatures considered here, but ex¬ 
hibit the logarithmically-diverging growth with decreas¬ 
ing temperature that is standard for a weak-coupling 
Cooper instability. As in the weak-coupling calculation, 
pairing between majority particles clearly dominates over 
minority pairing. 

Phase diagrams at densities n = 0.8 and 0.9 look very 
similar to the quarter filled case, indicating that the 
FFLO instability is not very sensitive to particle den¬ 
sity. For nearly half-filled bands, density-wave instabili¬ 
ties may however become relevant due to nesting in the 
particle-hole channel; the full phase diagram in the vicin¬ 
ity of half filling is therefore left for further studies. We 
have not systematically studied other interactions, but 
expect I//t = —4 to be close to the optimal case for ob¬ 
servation of FFLO order: At smaller |I/|, transition tem¬ 
peratures and the width of the FFLO instability will de¬ 
crease quickly im, whereas in the strong-coupling regime 
a (coherent or incoherent) mixture of tightly bound pairs 
and unpaired excess particles may be more stable than 
an FFLO phase. Note that a well-known particle-hole 
transformation relates the attractive and repulsive Hub¬ 
bard models to each other [miis]. Under this transfer- 
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mation, the magnetization nP assumes the role of doping 
X = n — 1 and vice versa, whereas the FFLO instability 
translates into an instability to a striped phase. 

Some questions concerning the extent and character of 
the FFLO state cannot be answered definitively by our 
study since we cannot enter the broken-symmetry phase. 
This concerns in particular the type of order (single-Q 
vs. multi-Q) and possible transitions between different 
ordered phases. We have not detected any hints of phase 
separation in the polarization vs. magnetic field curves, 
but we cannot conclusively rule this scenario out — even 
though previous studies on 2D lattices generally find di¬ 
rect and continuous transitions to the FFLO phase. 

In summary, we have presented the first well-controlled 
numerical evidence for the presence of a Fermi liquid 
instability towards FFLO order in the spin-imbalanced 
phase diagram of attractively interacting fermions on a 
2D lattice. For moderate on-site interaction Ujt = —4, 
the instability is present in an extended region of the 
temperature-polarization plane. The largest tempera¬ 
tures where this instability is observable are roughly by 
a factor of two smaller than the Kosterlitz-Thouless tran¬ 
sition temperature in the corresponding spin-balanced 
system, similar to DMFT results for anisotropic optical 
lattices [34]. At large polarization there does not seem 
to be any singlet superfluid order and triplet pairing is 
found at exponentially low temperatures. Our quantita¬ 
tive predictions may be checked by experiments with cold 
atoms in optical lattices, where the breaking of transla¬ 
tional symmetry can be observed by in-situ imaging [50]- 
[52] and the pair-momentum distribution by time-of-flight 
measurements [53H56] or in noise correlations [SilllH]. 
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